% twoD_extractor.m  2/12/2014  Parker MacCready
%
% this is for extracting the results of ROMS simulations
% it is focused on time series of 2D fields


clear;
% &&&&&&&&&&& USER EDIT THIS &&&&&&&&&&&&&&
doLocal = true;
if doLocal
    indir = '~/Documents/';
else
    indir = '/pmraid1/daves/runs/'; % for skua
end

whichCase = 'meanSurfaceVel';

%whichCase = 'max_bottom_vel';
% &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

addpath('../alpha')
[Tdir] = toolstart; % locate things and set paths

%% create output directory if needed
outdir = [Tdir.output,'twoD_out/'];
if ~exist(outdir,'dir'); mkdir(outdir); end;

%% choose run to work on and set the basename
disp('********* twoD_extractor.m *******************')
disp(' ')
if doLocal % interactive version
    disp('** Select folder where history files are **');
    pth = uigetdir(indir);
else % hardwired version
    pth = [indir,'salish_2006_4/OUT'];
end
ind = strfind(pth,'/');
basename = pth(ind(end-1)+1:ind(end)-1);
disp(' '); disp(['basename = ',basename]); disp(' ')

% get the run definition
run = roms_createRunDef('my_run',pth);

% get the year
year = datestr(run.his.nctime(1),'yyyy');
year = str2num(year);

% && USER EDIT &&
t0 = datenum(year,6,18);
t1 = datenum(year,6,24);
% &&&&&&&&&&
tag0 = num2str(round(t0));
tag1 = num2str(round(t1));

% prepare output arrays
lon = run.grid.lon;
lat = run.grid.lat;
mask = run.grid.mask;
lonu = run.grid.lonu;
latu = run.grid.latu;
lonv = run.grid.lonv;
latv = run.grid.latv;
[nr,nc] = size(lon);
max_speed = nan(nr,nc);
max_u = nan(nr,nc);
max_v = nan(nr,nc);
min_u = nan(nr,nc);
min_v = nan(nr,nc);
mean_u = zeros(nr,nc);
mean_v = zeros(nr,nc);

%%
switch whichCase
    case 'max_bottom_vel' % max and min near-bottom speed
        varlist = {'H','mab','max_u','max_v','max_speed', ...
            'lon','lat','min_u','min_v','mask'};
        for tt = run.his.ncn(1):run.his.ncn(end)
            tt_str = ['0000',num2str(tt)];
            tt_str = tt_str(end-3:end);
            % get the correct infile
            [infile] = [run.his.dirname,run.his.basename,tt_str,'.nc'];
            if tt == run.his.ncn(1)
                H = run.grid.H;
                z = roms_z(H,0*H,run.grid.cs);
                mab = squeeze(z(1,:,:)) + H; % mab of the deepest rho point
            end
            if mod(tt,5)==0; disp(['  infile = ',infile]); end;
            % near-bottom velocities
            u = nc_varget(infile,'u',[0 0 0 0],[1 1 -1 -1]);
            v = nc_varget(infile,'v',[0 0 0 0],[1 1 -1 -1]);
            ur = interp2(lonu,latu,u,lon,lat);
            vr = interp2(lonv,latv,v,lon,lat);
            speed = sqrt(ur.^2 + vr.^2);
            max_speed = max(max_speed,speed);
            max_u = max(max_u,ur);
            max_v = max(max_v,vr);
            min_u = min(min_u,ur);
            min_v = min(min_v,vr);
        end
        
    case 'meanSurfaceVel'
        varlist = {'H','mean_u','mean_v','lon','lat','mask'};
        it0 = dsearchn(run.his.nctime',t0);
        it1 = dsearchn(run.his.nctime',t1);
        ncn = run.his.ncn(it0:it1);
        NT = length(ncn);
        NZ = length(run.grid.cs);
        for tt = ncn(1):ncn(end)
            tt_str = ['0000',num2str(tt)];
            tt_str = tt_str(end-3:end);
            % get the correct infile
            [infile] = [run.his.dirname,run.his.basename,tt_str,'.nc'];
            if tt == ncn(1)
                H = run.grid.H;
            end
            if mod(tt,5)==0; disp(['  infile = ',infile]); end;
            u = nc_varget(infile,'u',[0 NZ-1 0 0],[1 1 -1 -1]);
            v = nc_varget(infile,'v',[0 NZ-1 0 0],[1 1 -1 -1]);
            ur = interp2(lonu,latu,u,lon,lat);
            vr = interp2(lonv,latv,v,lon,lat);
            mean_u = mean_u + ur/NT;
            mean_v = mean_v + vr/NT;
        end
end

%%  save the results as a netcdf file
ncfile_out = [outdir,whichCase,'_',tag0,'_',tag1,'.nc'];
disp(['Creating ' ncfile_out]);
my_mode = bitor ( nc_clobber_mode, nc_64bit_offset_mode );
nc_create_empty ( ncfile_out, my_mode );
% add dimensions
nx = length(lon(1,:));
ny = length(lat(:,1));
nc_add_dimension( ncfile_out, 'xi_rho', nx);
nc_add_dimension( ncfile_out, 'eta_rho', ny);
for ii = 1:length(varlist)
    varname = varlist{ii};
    varstruct.Name = varname;
    switch varname
        case 'H'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'bathymetry';
            units = 'meter';
            coordinates = 'lat lon';
            field = 'bath, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
        case 'mab'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'meters-above-bottom of speeds';
            units = 'meter';
            coordinates = 'lat lon';
            field = 'mab, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
            
        case 'mean_u'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'mean surface east-west speed';
            units = 'meters second-1';
            coordinates = 'lat lon';
            field = 'mean_u, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
        case 'mean_v'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'mean surface north-south speed';
            units = 'meters second-1';
            coordinates = 'lat lon';
            field = 'mean_u, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
            
        case 'max_u'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'max near-bottom east-west speed';
            units = 'meters second-1';
            coordinates = 'lat lon';
            field = 'max_u, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
        case 'max_v'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'max near-bottom north-south speed';
            units = 'meters second-1';
            coordinates = 'lat lon';
            field = 'max_u, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
        case 'min_u'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'min near-bottom east-west speed';
            units = 'meters second-1';
            coordinates = 'lat lon';
            field = 'min_u, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
        case 'min_v'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'min near-bottom north-south speed';
            units = 'meters second-1';
            coordinates = 'lat lon';
            field = 'min_u, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
        case 'max_speed'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'max near-bottom speed';
            units = 'meters second-1';
            coordinates = 'lat lon';
            field = 'max_speed, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
        case 'lon'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'longitude of RHO-points';
            units = 'degree_east';
            field = 'lon, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
        case 'lat'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'latitude';
            units = 'degree_north';
            field = 'lat, scalar';
            varstruct.Attribute = struct('Name', ...
                {'long_name','units','field'}, ...
                'Value',{long_name,units,field});
        case 'mask_rho'
            varstruct.Nctype = 'double';
            varstruct.Dimension = {'eta_rho';'xi_rho'};
            long_name = 'mask on RHO-points';
            option_0 = 'land';
            option_1 = 'water';
            coordinates = 'lat lon';
            varstruct.Attribute = struct('Name', ...
                {'long_name','option_0','option_1','coordinates'}, ...
                'Value',{long_name,option_0,option_1,coordinates});
    end
    nc_addvar(ncfile_out,varstruct);
    nc_varput(ncfile_out, varname, eval(varname));
end
